#!/usr/bin/python

#  This script takes two arguments, the input text integral file name and
#  the output text integral file name.
#  The scripts writes to the output file the hamiltonian integral
#  information for two non-interacting (infinitely separated) copies
#  of the system specified by the input file.

import sys
import string



def get_item(x, words, name_of_item):

  # raise an error of there are no more entries
  if len(words) == 0:
    raise RuntimeError("reached end of file before reading %s" % name_of_item)

  # integer
  if type(x) is int:
    return int(words.pop(0))

  # float
  if type(x) is float:
    return float(words.pop(0))

  # unsupported type
  raise TypeError("expected the type of x to be int or float")




if len(sys.argv) != 3:
  print "correct usage:    double_text_ints.py  input_file  output_file"
  sys.exit(1)

# read text from the input file
fin = open(sys.argv[1], "r")
words = []
for line in fin:
  words.extend(string.split(line))
fin.close()
i = len(words) - 1
while i >= 0:
  if len(words[i]) == 0:
    words.pop(i)
  i -= 1

# get number of orbitals
norb = get_item(int(1), words, "number of orbitals")
#print "read norb = %i" % norb

# get number of irreps
nirr = get_item(int(1), words, "number of irreps")
#print "read nirr = %i" % nirr

# get orbital irreps
orb_irrep = [ int(-1) for i in range(norb) ]
for i in range(norb):
  orb_irrep[i] = get_item(int(1), words, "orbital irrep for orbital %i" % i)
#print "read orbital irreps ", orb_irrep

# get core/nuclear energy
e_cor_nuc = get_item(float(1.0), words, "core/nuclear energy")
#print "read core/nuc energy %20.12f" % e_cor_nuc

# overlap integrals
ovi = [ [ float(0.0) for i in range(norb) ] for j in range(norb) ]
for p in range(norb*norb):
  i = get_item(int(1), words, "1st ovi index")
  j = get_item(int(1), words, "2nd ovi index")
  ovi[i][j] = get_item(float(1.0), words, "ovi[%i][%i]" %(i,j))

# one electron integrals
oei = [ [ float(0.0) for i in range(norb) ] for j in range(norb) ]
for p in range(norb*norb):
  i = get_item(int(1), words, "1st oei index")
  j = get_item(int(1), words, "2nd oei index")
  oei[i][j] = get_item(float(1.0), words, "oei[%i][%i]" %(i,j))

# two electron integrals
nread = 0
tei = []
try:
  while True:
    nread = 0
    i = get_item(int(1), words, "1st tei index")
    nread += 1
    j = get_item(int(1), words, "2nd tei index")
    nread += 1
    k = get_item(int(1), words, "3rd tei index")
    nread += 1
    l = get_item(int(1), words, "4th tei index")
    nread += 1
    val = get_item(float(1.0), words, "4th tei index")
    nread += 1
    tei.append([i,j,k,l,val])
except RuntimeError:
  if nread != 0:
    raise RuntimeError("reached end of file while only partially finished reading a two electron integral")

# open output file
fout = open(sys.argv[2], "w")

# write number of orbitals
fout.write("%i\n\n" % (2*norb))

# write number of irreps
fout.write("%i\n\n" % (nirr))

# write orbitals irreps
for i in orb_irrep:
  fout.write(" %i" % (i))
for i in orb_irrep:
  fout.write(" %i" % (i))
fout.write("\n\n")

# write core/nuc energy
fout.write("%20.12f\n\n" % (2.0*e_cor_nuc))

# write overlap integrals
for i in range(2*norb):
  for j in range(2*norb):
    if i < norb and j < norb:
      fout.write("%4i  %4i  %20.12f\n" % ( i, j, ovi[i][j] ) )
    elif i >= norb and j >= norb:
      fout.write("%4i  %4i  %20.12f\n" % ( i, j, ovi[i-norb][j-norb] ) )
    else:
      fout.write("%4i  %4i  %20.12f\n" % ( i, j, 0.0 ) )
fout.write("\n")

# write one electron integrals
for i in range(2*norb):
  for j in range(2*norb):
    if i < norb and j < norb:
      fout.write("%4i  %4i  %20.12f\n" % ( i, j, oei[i][j] ) )
    elif i >= norb and j >= norb:
      fout.write("%4i  %4i  %20.12f\n" % ( i, j, oei[i-norb][j-norb] ) )
    else:
      fout.write("%4i  %4i  %20.12f\n" % ( i, j, 0.0 ) )
fout.write("\n")

# write two electron integrals
for t in tei:
  fout.write("%4i  %4i  %4i  %4i  %20.12f\n" % ( t[0], t[1], t[2], t[3], t[4] ) )
for t in tei:
  fout.write("%4i  %4i  %4i  %4i  %20.12f\n" % ( t[0]+norb, t[1]+norb, t[2]+norb, t[3]+norb, t[4] ) )
fout.write("\n")

# close the output file
fout.close()
